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Abstract 


A  modified  form  of  the  MODFLOW-2000  code  provides  a  tool  to  estimate  a  hydraulic  conductivity  distribution 
near  a  borehole.  The  modified  code  allows  the  effects  of  wellbore  screen  intervals,  screen  clogging,  and  disturbed 
zone  skin  to  be  taken  into  account.  The  primary  application  has  been  in  an  investigation  of  electromagnetic 
borehole  flow  meter  measurement  interpretation. 

MODFLOW-2000  was  modified  to  simulate  a  two-dimensional  cylindrical  geometry.  The  two  dimensions  axe 
vertical  and  radial  with  angular  symmetry.  The  geometry  is  similar  to  the  RADMOD  preprocessor  to  MOD- 
FLOW,  but  is  more  flexible  in  that  the  grid  spacing  in  the  radial  direction  may  be  non-regular  and  the  hydraulic 
conductivity  of  a  layer  need  not  be  constant.  This  modification  includes  sensitivity  calculations  so  that  parameter 
estimation  can  be  performed. 

To  support  parameter  estimation  using  flow  meters,  a  new  measurement  category  of  internal  flow  observations 
was  introduced.  Internal  flow  observations  occur  at  locations  that  are  fully  contained  within  the  model  domain. 
The  standard  MODFLOW-2000  code  is  limited  to  flow  observations  across  the  model  domain  boundaries.  The  new 
internal  flow  measurement  observations  are  compatible  with  the  normal  MODFLOW-2000  Cartesian  geometry 
and  with  the  cylindrical  geometry  modifications. 

A  third  modification  moves  the  node  location  of  a  cell  containing  the  water  table  from  the  center  of  the  cell 
to  the  location  of  the  water  table.  The  vertical  position  of  the  node  dynamically  follows  changes  in  the  water 
table  position.  This  change  improves  the  simulation  of  water  table  dynamics  for  both  cylindrical  geometry  and 
Cartesian  geometry  simulations. 

Together  these  modifications  provide  a  tool  for  modeling  pumping  from  a  wellbore  with  multiple  screened 
zones  in  a  layered  aquifer  with  an  azimuthally  uniform  disturbed  zone.  If  the  wellbore  is  included,  the  pumping 
rate  can  be  applied  to  the  upper  node  of  the  wellbore.  The  flow  into  or  from  the  formation  is  correctly  apportioned 
in  the  model.  Thus,  the  modeler  does  not  have  to  apportion  flow  into  individual  layers,  a  difficult  procedure  in 
the  case  of  a  disturbed  zone  that  results  in  a  positive  skin  effect. 
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1  Introduction 


The  finite-difference  formulation  for  cylindrical  geometry  used  by  the  RADMOD  code  (Reilly  and  Harbaugh,  1993) 
provides  a  capability  to  simulate  radial  flow  to  a  well.  However,  RADMOD  lacks  the  flexibility  to  simulate  radial 
variation  in  hydraulic  conductivity  caused  by  clogged  well  screens  or  formation  disturbances  near  the  borehole. 
We  need  this  flexibility  to  an  investigate  the  influence  of  such  heterogeneity  on  wellbore  flow  profiles  created  from 
electromagnetic  borehole  flow  meter  measurements. 

Briefly,  an  electromagnetic  borehole  flow  meter  measures  the  voltage  drop  across  an  annulus,  created  by  ions 
moving  through  a  magnetic  field  imposed  in  the  annulus  by  the  flowmeter.  The  voltage  drop  is  used  to  determine  the 
axial  flow  velocity  within  a  borehole  (Young  and  Pearson,  1995).  The  net  inflow  to  the  well  over  an  interval  can  be 
estimated  by  calculating  the  difference  in  vertical  flow  measured  at  the  two  ends  of  the  interval.  Prom  the  net  inflow 
one  can  deduce  estimates  of  the  hydraulic  conductivity  distribution  of  the  formation  near  the  borehole.  Common 
practice  at  this  time  is  to  assume  a  perfectly  layered  formation  resulting  in  conductivity  estimates  proportional  to 
the  ratio  of  inflow  rates  to  the  borehole  (Molz  et  al.,  1994).  At  the  Boise  Hydrogeophysical  Research  Site,  BHRS 
(Barrash  et  al.,  1999),  we  believe  that  the  wellbore  screens  may  be  partially  clogged  with  sand  resulting  in  a  lower 
hydraulic  conductivity  through  the  screen  compared  to  the  formation.  Our  initial  investigations  using  a  thermal 
simulator,  QUICKFIELD  (Tera  Analysis,  1999),  revealed  that  clogging  causes  a  distortion  of  the  flow  into  a  wellbore 
compared  to  the  flow  in  a  layered  formation,  potentially  biasing  the  interpretation  of  electromagnetic  flowmeter 
measurements. 

We  decided  to  replace  QUICKFIELD  in  our  analysis  because  a  lack  of  access  to  the  source  code  makes  it 
cumbersome  to  integrate  into  an  analysis  scheme  and  because  we  felt  a  fluid  flow  simulator  would  be  more  readily 
accepted  by  the  hydrologic  community.  Of  the  readily  available  non-commercial  software  we  examined,  RADMOD, 
and  the  similar  capabilities  of  RADFLOW  (Johnson  et  al.,  2002),  came  the  closest  to  fulfilling  our  needs.  However, 
we  could  not  specify  the  hydraulic  conductivity  of  the  screen  to  be  different  from  that  of  the  formation.  We  would 
also  have  to  base  the  geometry  of  the  formation  on  the  geometry  of  the  screen,  which  would  have  led  to  a  large  model. 
We  chose  to  modify  MODFLOW-2000  (Harbaugh  et  al.,  2000;  HiU  et  al.,  2000)  rather  than  RADMOD  because  of 
usefulness  to  us  of  the  inversion  capability  of  MODFLOW-2000.  This  report  is  limited  to  the  modifications  that 
were  made  to  MODFLOW-2000.  The  investigation  of  electromagnetic  borehole  flow  meter  interpretation  will  be 
documented  separately. 

With  our  modified  version  we  can  include  the  open  wellbore  explicitly  in  our  simulations  of  electromagnetic  flow 
meter  measurements.  In  order  to  provide  simulated  measurements  within  the  borehole,  a  new  observation  type 
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was  needed.  The  2001  version  of  MODFLOW-2000  can  be  used  to  simulate  observations  of  flow  across  the  model 
domain  boundaxies,  but  not  internal  flows  that  do  not  cross  a  boundary.  So,  we  added  an  internal  flow  observation 
type.  The  internal  flow  observation  capability  is  fully  integrated  with  the  other  observation  types  and  could  be  used 
independently  of  the  other  modifications  described  in  this  report. 

While  benchmarking  the  transient  behavior  of  the  cylindrical  geometry  modifications,  the  poor  response  of  the 
model  near  the  water  table  was  noted.  RADMOD  provides  a  better  simulation  of  water  table  dynamics  than  we 
achieved  with  our  initial  modifications.  Examining  the  RADMOD  code  revealed  that  the  outer  nodes  are  located  at 
the  domain  boundaries  in  a  RADMOD  model.  Having  nodes  on  the  model  edge  is  incompatible  with  MODFLOW- 
2000  but  is  compatible  with  earlier  versions  of  MODFLOW  using  the  General  Finite-Difference  Package.  General 
Finite-Difference  Package  allows  more  flexible  model  geometry  compared  to  standard  MODFLOW.  RADMOD  is  a 
preprocessor  that  supplies  an  input  file  for  the  General  Finite-Difference  Package.  As  we  show  below,  dynamically 
positioning  the  node  of  a  finite-difference  cell  at  the  water  table  leads  to  a  more  accurate  simulation  of  water  table 
dynamics. 

Since  our  objective  was  to  create  a  tool  for  our  investigation  of  borehole  flow  meter  measurements,  we  have  not 
attempted  to  make  the  cylindrical  geometry  model  compatible  with  all  of  the  MODFLOW-2000  packages.  While 
we  hope  to  completely  integrate  the  cylindrical  geometry  model  into  MODFLOW-2000,  we  will  not  do  so  unless  the 
USGS  expresses  an  interest.  The  current  limitations  of  our  modifications  are  listed  in  Section  5  of  this  report. 

The  exact  coding  changes  made  to  MODFLOW-2000  do  not  belong  in  this  report.  However,  the  changes  may  be  of 
interest  to  its  readers.  We  provide  Internet  access  to  these  modifications  through  http://cgiss.boisestate.edu/cgiss_pub.hi 
The  file  is  accessible  through  the  listing  of  this  technical  report,  CGISS  02-01.  The  file  is  a  compressed  tape  archive 
(tar)  file.  The  compression  was  performed  using  Gnu’s  gzip  utility.  The  file  contains  all  the  modified,  and  unmodified, 
files  needed  to  compile  the  code.  A  file  for  the  make  utility  is  also  supplied.  UNIX  diff  command  output  files,  in 
directory  diff,  contain  just  the  changes  made  to  Version  1.7  12/04/2001  of  MODFLOW-2000.  In  addition,  the  tar 
file  includes  the  two  test  cases  described  in  Sections  2.3  and  4.4.1.  Finally,  this  report  is  included.  The  file  contents 
may  be  updated  in  the  future  to  reflect  corrections  and  additional  compatibility  improvements. 
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2  Cylindrical  Geometry 

2.1  Formulation 

2.1.1  Conductance 

The  description  of  modifications  made  to  support  cylindrical  geometry  follows  the  development  of  the  RADMOD 
code  by  Reilly  and  Harbaugh  (1993).  The  modification  differs  from  RADMOD  with  the  relaxation  of  the  use  of  a 
constant  grid-expansion  factor.  In  this  section,  we  assume  the  reader  has  some  familiarity  with  MODFLOW,  but  not 
necessarily  with  RADMOD.  Figure  1  depicts  the  finite-difference  discretization  in  cylindrical  geometry.  The  model 
domain  is  composed  of  a  series  of  annular  rings  about  a  common  central  location.  The  modified  code  accommodates 
constant  layer  thicknesses,  A Z,  and  arbitrary  radial  spacing,  A R.  (Variable  layer  thickness  might  also  be  simulated 
correctly,  but  has  not  been  investigated.)  In  the  figure,  the  lines  separate  finite-difference  cells.  For  a  few  cells,  the 
cell  computational  nodes  are  shown  as  dots.  The  location  of  a  cell  boundary  is  specified  by  the  radial  coordinate  sr. 
The  radial  position  of  a  node  location  is  designated  r,  as  shown  for  the  last  column  in  the  figure.  The  inner  domain 
boundary  of  the  model  may  be  the  central  axis  (sri=0.0)  or  at  some  radial  distance  away  from  the  axis. 


z 


Figure  1:  Finite-difference  discretization  for  cylindrical  geometry. 

Each  cell  in  a  model  is  assigned  a  set  of  constant  properties.  The  MODFLOW  codes  transform  the  cell’s  assigned 
hydraulic  conductivity  or  transmissivity  into  hydraulic  conductance  terms  between  nodes.  Hydraulic  conductance 
includes  the  geometry  of  the  finite-difference  cells  in  the  determination  of  the  effective  ability  to  transmit  flow  between 
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nodes.  The  basic  equation  for  conductance  is 


(I) 


where:  K  is  hydraulic  conductivity,  A  is  the  area  of  the  common  boundary  between  cells  and  L  is  the  distance 
between  nodes.  Reilly  and  Harbaugh  present  the  hydraulic  conductance  in  the  radial  direction  between  two  nodes 

in  cylindrical  geometry  as 


Cr  = 


KT2nAz 


(2) 


where:  Kr  is  the  hydraulic  conductivity  between  the  nodes  located  at  radii  of  n  and  ri+iand  Az  is  the  axial  thickness 
of  the  layer.  The  nodes  are  the  “centers”  of  two  adjoining  cells.  The  cell  “center”  is  the  radial  location  that  separates 
the  cell  into  two  halves  of  equal  cross  sectional  area  normal  to  the  vertical  direction,  rather  than  the  location  with 
equal  distance  from  the  edges  of  the  cell.  Equation  2  is  based  on  constant  hydraulic  conductivity  within  each  layer. 

In  our  case,  we  have  constant  hydraulic  conductivity  within  each  cell,  but  each  cell  may  have  a  different  hydraulic 
conductivity.  Therefore,  we  calculate  conductance  for  each  half  cell  between  the  center  node  and  the  cell  boundary 
for  adjacent  cells  and  then  add  the  two  half  cell  conductances  in  series.  For  the  two  half  cell  conductances  we  have 


_  Kr.  2nAz 

~ '» (if) 


(3) 


and 


^n+i  — 


to(^) 


where,  sri  is  the  boundary  between  cell  i  and  cell  *  + 1.  The  conductance  between  the  nodes  becomes 

Cu  Crt+i 


Cr  = 


CTi  +  C, 


(4) 


n+i 


Note  that  if  Ku  =  Ku+1  =  Kr  then  (4)  is  equivalent  to  (2). 

Continuing  to  follow  Reilly  and  Harbaugh,  the  cell  “centers”  are  located  such  that  the  cross-sectional  areas  normal 
to  the  z  axis  of  the  two  half  cells  are  equal.  The  two  half  cells  are:  1)  less  than  the  center  position,  and  2)  greater 
than  the  center  position. 

i r  (rf  -  sr?)  =  n  ( srf+1  -  r?)  (5) 


Thus, 


U  =  \j\  (srf+  sr?+1) 


(6) 
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Whereas  RADMOD  uses  a  constant  expansion  ratio  to  define  sr,  we  use  the  standard  MODFLOW  spacing  along 
rows,  DELR,  to  define  the  cell  boundaries.  One  absolute  position,  the  radius  of  the  inner  boundary  of  the  first 
column,  is  also  required  to  fully  define  the  sr  values.  The  definition  of  the  vertical  conductance  follows  directly  from 
(1)  with  the  definition  of  cross-sectional  area  as 

Ai  =  7T  (sr?  -  sr?+1 )  (7) 


2.1.2  Conductance  and  the  Sensitivity  Process 

The  inversion  capability  of  MODFLOW-2000  requires  the  determination  of  the  sensitivity  of  observations  to  pa.ra.m- 
eters.  By  the  chain  rule  of  calculus  the  sensitivity  of  observations  to  parameters  can  be  calculated  as  the  sensitivity 
of  the  observations  to  conductance  multiplied  by  the  sensitivity  of  the  conductance  to  the  parameters.  The  required 
modification  of  the  MODFLOW-2000  sensitivity  process  is  quite  limited  because  the  sensitivities  of  the  observations 
to  conductance  remain  unchanged. 

As  an  example,  consider  an  internal  flow  observation  that  is  defined  as  the  radial  flow  from  cell  i  to  *  +  1.  The 
observation  can  be  written  as 

Qo  =  Cr  (hj+i  -  hi)  (8) 

The  sensitivity  to  a  hydraulic  conductivity  parameter,  6,  is 


dQo 

db 


dCr 

db 


(hi+i  —  hi)  4-  Cr 


(  dhj+i 

V  9Cr 


dhi 

dCr 


\  dCr 

)  db 


(9) 


The  terms  ^  do  not  require  modification  for  cylindrical  geometry.  We  only  need  to  provide  a  calculation  of 
where  x  represents  either  the  vertical  or  radial  direction.  By  the  chain  rule 


dCr  <%+1  dCri  C%  dCri+1 

db  (Cri+Cri+1)2  db  (Cri+Cri+1)2  db 


dCTi  _  2ttAz  9Cri+1  _  2-kAz 

db  "  in(£*)  ,an  56  -  Jajfi 

Likewise, 

dCVi  Ai 

db  |  (A  zu  +  A  zi) 

where,  u  and  l  refer  to  the  thickness  of  the  upper  layer  and  lower  layers  respectively. 


(10) 


(11) 


(12) 
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2.1.3  Storage  Terms 

The  storage  terms  are  simply  volume,  Vi,  times  specific  storage,  Ss, 

SSi  =  S,Vi 
=  SsAiAz 


or  area,  A,,  times  specific  yield,  5„,  in  the  case  of  imconfined  layers. 

SYi  =  SyAi  O'1) 

where  At  is  defined  by  (7).  The  sensitivities  of  storage  terms  to  specific  storage  parameters  or  specific  yield  parameters 
are  the  volume  or  area  respectively. 

2.1.4  Horizontal  Flow  Barriers 

The  horizontal  flow  barrier  package  allows  a  thin  low  permeability  zone  to  be  incorporated  into  the  model  without 
having  to  explicitly  represent  the  geometry  of  the  zone.  For  a  borehole  flowmeter  investigation,  the  horizontal  flow 
barriers  can  be  used  to  simulate  well  casing  and  to  represent  clogging  of  the  wellbore  screen.  Within  MODFLOW- 
2000,  the  barrier  is  treated  as  a  conductance  at  the  boundary  between  cells  that  is  added  in  series  to  the  horizontal 
conductance.  The  data  input  for  a  horizontal  barrier  is  in  the  form  of  hydraulic  conductivity  divided  by  the  barrier 
thickness,  K/L,  hence  the  thickness  need  not  be  explicitly  defined.  The  conductance  of  the  barrier,  Cb,  is  the 
cross-sectional  area  normal  to  the  flow  direction  multiplied  by  the  data  input. 

Cb  =  Ab-  (K/L)  (15) 

In  the  radial  flow  case,  the  cross-sectional  is  Ab  =  2tt  srj  A z,  replacing  DELR-Az  or  DELC-Az  for  both  the  flow  and 
sensitivity  processes. 

2.2  Implementation 

The  cylindrical  geometry  option  is  specified  by  setting  the  conductance-averaging  flag,  LAYAVG,  in  the  LPF  file  to 
3  for  all  layers.  Setting  LAYAVG  to  3  for  some  but  not  all  of  the  layers  will  generate  an  error  message.  In  addition, 
if  the  cylindrical  geometry  option  is  used,  the  option  XSECTION  must  be  set  in  the  BAS6  file.  The  radial  distance 
of  the  innermost  boundary  is  read  from  a  file  type  CGEO,  IUNIT(50),  as  specified  in  the  name  file.  (See  Appendix 
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A.l).  The  radial  distance  of  the  innermost  boundary  is  the  only  datum  in  this  file.  Consistent  with  the  XSECTION 
specification,  the  DELR  array  in  the  DIS  file  is  used  to  specify  radial  cell  widths.  In  the  DIS  file,  the  number  of  rows, 
NROW,  must  be  set  to  1  and  the  width  of  the  columns,  DELC,  must  be  set  to  1.0.  No  automated  error  checking  of 
NCOL  or  DELC  is  provided. 

Four  internal  variables  are  calculated  within  MODFLOW-2000  to  implement  the  cylindrical  geometry  option. 
These  me  the  radial  conductance,  CR,  the  vertical  conductance,  CV,  the  cell  hydraulic  capacity,  SCI,  and,  for 
convertible  cells,  the  cell  water  table  yield,  SC2.  The  calculations  of  all  of  these  variables  are  controlled  by  a  new 
subroutine  called  SGWF1LPF1RADIAL  that  has  been  added  to  the  gwfllpfl.f  file.  The  cross-sectional  areas  and 
cell  node  radial  positions  are  also  calculated  and  stored  within  SGWF1LPF1RADIAL  using  the  ALLOCATEABLE 
variable  attribute.  Thus,  compilation  of  the  modified  code  is  limited  to  FORTRAN  90  and  95  compilers.  The  advan¬ 
tage  of  isolating  these  calculations  and  variables  within  the  SGWF1LPF1RADIAL  subroutine  is  that  perturbation 
of  the  standard  MODFLOW-2000  code  is  limited. 

A  second  subroutine,  SSEN1LPF1RADIAL,  is  used  to  calculate  the  derivatives  of  CR,  CV,  SCI,  and  SC2  with 
respect  to  parameters.  This  subroutine  is  also  called  by  the  horizontal  flow  barrier  package  to  obtain  cross-sectional 
areas  normal  to  the  radial  direction.  As  in  the  SGWF1LPF1RADIAL  subroutine,  the  cross-sectional  areas  normal  to 
the  vertical  direction  and  the  cell  node  positions  are  calculated  and  stored.  This  small  duplication  of  storage  space 
is  inefficient  but  limits  the  perturbation  of  the  original  code. 

2.3  Testing 

We  follow  Reilly  and  Harbaugh’s  (1993)  lead  by  using  the  test  case  presented  in  the  RADMOD  documentation:  flow 
to  a  partially  penetrating  well  in  a  homogeneous  aquifer  of  infinite  extent  with  anisotropic  hydraulic  conductivity. 
The  WTAQ3  (Barlow  and  Moench,  1999)  numerical  implementation  of  an  analytic  solution  (Moench,  1997)  for  a 
homogeneous  water-table  aquifer  response  to  pumping  from  a  partially  penetrating  finite-diameter  well  was  used  to 
define  the  reference  drawdown.  Moench’s  analysis  is  an  extension  of  the  analytical  solution  (Neuman,  1974)  that. 
was  used  to  benchmark  the  RADMOD  preprocessor.  The  reference  case  includes  pumping  from  a  25  foot  interval 
at  the  bottom  of  a  100  foot  thick  homogeneous  aquifer  of  infinite  extent.  The  aquifer  has  a  horizontal  hydraulic 
conductivity  of  100  ft/day  and  a  vertical  hydraulic  conductivity  of  10  ft/day.  The  transmissivity,  T,  of  the  aquifer 
is  10,000  ft2 /day.  The  specific  storage  is  5.xl0-6  per  ft  and  the  specific  yield  is  0.2.  A  pumping  rate,  Q,  of  125,670 
ft8 /day  is  used,  which  makes  the  dimensionless  drawdown  equivalent  to  the  drawdown  in  units  of  feet.  The  well  has 
a  radius  of  0.936  ft. 
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Figure  2:  Schematic  of  test  model  for  cylindrical  geometry. 

Note  that  the  horizontal  axis  is  a  logarithmic  scale.  The  dots  represent 
node  positions.  They  are  spaced  evenly  along  the  log  axis.  After  (Reilly  and  Harbaugh,  1993). 


As  shown  in  Figure  2,  the  test  case  has  11  layers  and  40  columns.  The  radial  position  of  the  nodes  increases 
exponentially  with  a  constant  growth  multiple  of  1.5.  In  Figure  2,  the  layer  numbers  are  indicated  on  the  vertical  axis 
and  the  column  numbers  are  shown  above  the  drawing.  Note  that  the  horizontal  axis  is  drawn  with  a  logarithmic 
scale  making  the  nodes  appear  to  have  equal  separation.  The  nodes  used  to  simulate  pumping  withdrawal  are  the 
three  lower  nodes  of  Column  1,  shown  in  the  lower  left  within  a  box.  The  pumping  rate  is  apportioned  between  these 
nodes  by  the  thickness  of  the  layers.  Three  nodes,  at  locations  Layer  6  Column  8,  Layer  7  Column  10,  and  Layer  3 
Column  2,  are  circled.  The  calculated  drawdown  at  these  nodes  is  used  to  compare  the  cylindrical  geometry  results 
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with  the  WTAQ  solution. 

To  prepare  for  the  test  case,  RADMOD  was  modified  to  report  the  radial  separation  distance  between  cell 
boundaries.  RADMOD  output  was  then  used  as  the  DELR  input  data.  Changes  were  made  to  the  RADMOD  model 
in  that  the  node  closest  to  the  axial  center  is  the  center  of  a  cell.  Thus,  the  innermost  boundary  of  the  MODFLOW- 
2000  model  is  located  at  a  radial  distance  of  0.708  rather  than  0.936,  which  positions  the  nodes  of  column  1  at  the 
same  position  in  both  models.  Likewise,  the  last  nodes  of  the  last  column  are  positioned  at  the  center  of  the  column, 
thus  extending  the  outer  boundary  of  the  model.  An  inconsequential  impact  of  these  differences  is  that  the  vertical 
cross-sectional  area  for  flow  is  larger  for  the  first  and  last  columns,  resulting  in  reduced  drawdown  near  the  well 
during  the  first  second  of  pumping. 

We  have  retained  the  positioning  of  the  nodes  at  the  top  layer  along  the  top  boundary,  at  least  at  the  beginning 
of  the  pumping  period.  The  MODFLOW-2000  code  has  been  modified  so  that  the  location  of  the  node  in  a  cell  that 
contains  a  water  table  is  dynamically  positioned  at  the  water  table.  Section  4  describes  in  detail  the  rationale  for 
positioning  of  the  upper  nodes  at  the  water  table.  The  change  has  a  significant  impact  on  the  model  performance. 
The  nodes  in  Layer  11,  the  bottom  layer,  are  located  at  the  center  of  the  layer  in  the  MODFLOW-2000  model  rather 
than  the  lower  boundary  as  in  the  RADMOD  model.  In  the  MODFLOW-2000  model,  the  nodes  of  the  bottom  layer 
are  closer  to  the  nodes  of  the  neighboring  layer  by  2.5  ft  compared  to  the  separation  distance  of  the  RADMOD 
model.  As  with  the  changes  to  Columns  1  and  40,  we  believe  that  this  change  has  little  effect  on  the  simulation. 

Figure  3  presents  the  comparison  of  dimensionless  drawdown  calculated  for  the  three  highlighted  nodes  of  Figure 
2.  The  MODFLOW-2000  calculations  are  shown  as  dashed  lines,  with  the  WTAQ  calculations  as  solid  lines.  The 
dimensionless  drawdown  is  calculated  from  dimensional  drawdown  as 


_  4t xT 

d*=-q‘ 


where  s  is  the  drawdown  in  units  of  feet.  The  abscissa  is  dimensionless  time  -  calculated  as 


T  Tt 
~  0  -2 


Str{ 


6,8 


(16) 


(17) 


where:  t  is  the  time  in  days,  St  is  storativity,  5xl0"4,  and  re, 8  is  16  ft,  the  radial  position  of  the  node  located  at 
Layer  6  Column  8.  The  comparisons  reveal  a  slight  under-prediction  of  drawdown  by  MODFLOW-2000. 

The  modified  MODFLOW-2000  code  allows  wellbore  dynamics  to  be  simulated,  which  is  not  possible  using  the 
RADMOD  preprocessor.  The  test  case  was  changed  so  that  the  first  column  of  the  model  represents  the  wellbore. 
For  this  new  case,  the  location  of  the  inner  boundary  is  the  z  axis  (sri=0.0  ft)  and  the  location  of  the  boundary 
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between  the  first  and  second  column  is  0.936  ft.  The  hydraulic  conductivity  of  the  wellbore,  Column  1,  was  set  to  an 
arbitrarily  large  value  of  10“ft/day.  The  specific  storage  of  the  cells  hr  column  one  was  reduced  to  6.0*10-“  per  ft 
and  the  specific  yield  was  increased  to  1.0.  The  homontal  flow  barrier  package  was  used  to  simulate  the  impermeable 
borehole  casing  by  introducing  a  10  “  per  day  conductance  barrier  between  the  firs,  and  second  column  of  layers  , 
through  8.  The  WEL  file  was  changed  to  specify  pumping  of  the  entire  125,670  ft“/day  from  Layer  5,  Column  1. 
Layer  5  was  ch«en  because  a  drawdown  of  40  ft  occurs  in  the  well,  cauang  the  upper  three  layers  of  the  wellbore  to 
go  dry.  In  the  WTAQ  input,  the  simulation  was  changed  from  an  infinitesimal  well  to  a  well  with  an  0.936  ft  innde 

radius. 

Figure  4  presents  the  dimensionless  curves  for  each  of  the  three  comparison  locations.  The  WTAQ  calculations 
for  an  infimtesimal  radius  pumping  well  are  drawn  as  dotted  fines.  The  solid  fines  are  for  the  finite  dimensional 
wellbore  and  the  dashed  lines  are  from  the  MODFLOW-2600  simulations. 


Figure  3:  Comparison  of  cylindrical  geometry  model  (dashed)  ft)  WTAQ  calculations  (solid). 


Ficure  4:  Comparison  of  drawdown  with  wellbore  storage. 

The  solid  lines  are  WTAQ  results,  dashed  lines  are  MODFLOW-2000, 
and  the  dotted  lines  are  from  WTAQ  without  wellbore  storage 
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3  Internal  Flow  Observations 


3.1  Introduction 

Internal  flow  observations  are  calculations  of  the  flow  between  cells  of  the  model.  Groups  of  cells  may  be  defined 
to  specify  flow  across  long  cross-sections  through  a  model.  This  could  be  used  as  an  alternative  to  calculating  flow 
through  portions  of  a  model  from  the  cell-by-cell  flow  file.  The  internal  flow  observation  calculations  axe  performed 
in  double  precision  so  the  internal  flow  observations  may  be  more  accurate  than  post-processing  of  the  MODFLOW 
results.  Our  use  is  in  conjunction  with  the  cylindrical  geometry  model  to  simulate  borehole  flow  meter  measurements. 

The  code  modification  for  the  internal  flow  observations  is  distinct  from  the  cylindrical  geometry  modification. 
Internal  flow  observations  do  not  require  the  use  of  cylindrical  geometry  nor  are  the  internal  flow  observations  required 
to  use  cylindrical  geometry. 


3.2  Formulation 

Internal  flow  observations  are  the  sum  of  flows  crossing  specific  cell  boundaries.  The  flow  is  calculated  as 


Qi  —  Cxi  *  (hi  “  hj) 


(18) 


where  CXi  is  the  conductance  in  the  row,  column,  radial,  or  vertical  direction  between  contiguous  cells  i  and  j,  and 
hi  is  the  hydraulic  head  in  the  cell. 

The  sensitivity  of  the  internal  flow  observations  to  parameter  changes  is 


dQi 

db 


dCxi 

8b 


{hi  —  ft»+x)  +  CXi 


(dhi 

V  db 


dhj+i  \ 
db  ) 


(19) 


Equation  9  described  the  same  sensitivity  in  terms  of 


3.3  Implementation 

Internal  flow  observations  are  activated  by  supplying  a  data  file  designated  as  type  FLOB  in  the  name  file.  The 
format  of  the  data  follows  that  of  the  general  head  boundary  observations  but  is  not  identical.  Internally  within  the 
code,  the  existence  of  internal  flow  observations  is  recognized  by  the  non-zero  status  of  variable  IUNIT(42).  A  set 
of  subroutines  following  the  structure  of  the  general  head  boundary  observations  is  supplied  in  file  obslflwl.f.  The 
subroutine  naming  convention  is  followed,  with  FLW  substituted  for  GHB.  The  data  storage  of  the  internal  flow 
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observations  has  been  integrated  with  the  other  flow  observations.  It  is  assumed  that  the  internal  flow  observations 
are  the  last  type  of  flow  obsecrations  added  to  the  list  of  observations.  This  may  need  to  be  changed  if  additional 

flow  observation  types  are  added  to  the  code. 

Formal  input  instructions  for  the  internal  flow  observations  axe  included  in  Appendix  A.2.  The  input  differs 
slightly  from  the  general  head  boundary  observation  input  because  the  internal  flow  observations  require  the  desig¬ 
nation  of  two  cells  rather  than  a  single  cell.  Item  5  of  the  input  has  been  changed  so  that  seven  data  items  axe  read 
for  each  component  of  an  observation  group.  The  layer,  row,  and  column  of  the  receiving  cell  axe  given  first.  The 
direction  of  flow  is  defined  so  that  flow  toward  the  receiving  cell  is  positive.  Next,  the  layer,  row  and  column  of  the 
discharging  cell  are  listed.  Positive  flow  is  away  from  the  discharging  cell.  The  seventh  data  item  read  is  the  weight 
factor.  The  weight  factor  is  a  multiplying  factor  to  allow  a  fraction  of  the  flow  between  individual  cells  to  contribute 
to  the  flow  measurement.  The  weight  factor  provides  a  correction  factor  if  the  location  of  the  flow  measurement  does 
not  coincide  perfectly  with  the  cell  geometry.  For  example,  if  the  grid  size  is  larger  than  the  measurement  interval, 
then  a  fraction  of  the  flow  from  cell  to  cell  should  be  calculated.  If  the  desired  flow  occurred  along  some  interior 
line  of  a  cell,  then  a  weighted  average  of  the  flow  across  the  cells  boundaries  might  be  used.  Multiple  cell  pairs  may 
be  specified  for  an  individual  flow  observation.  The  consistency  of  each  cell  pair  is  checked  to  ensure  that  they  are 
adjacent.  No  automated  check  of  the  continuity  of  multiple  cells  in  an  observation  is  performed. 

3.4  Example:  Borehole  Flow  Meter  Observation 

figure  5  depicts  the  region  near  the  pumped  eectiou  of  the  second  test  case  of  section  2.3.  The  figure  shows  the 
innermost  four  columns  and  bottom  five  layers.  In  this  test  case,  the  borehole  is  explicitly  included  in  the  simulation 
as  shown  by  the  cross-hatching  of  Column  1.  A  thick  line  separating  Columns  1  and  2  of  Layers  7  and  8  depicts  the 
casing  simulated  by  an  effectively  impermeable  horizontal  flow  barrier.  The  three  circled  regions  highlight  the  cell 
boundaries  we  now  use  to  represent  borehole  flow  meter  measurement  locations. 

In  the  FLOB  file,  we  specify  18  internal  flow  observations,  composed  of  one  cell  pair  for  each  of  the  three  locations. 
Observations  are  requested  at  five  times,  one  per  decade  of  dimensionless  time  from  1.0  (1.1  s)  to  104  (0.13  days). 
The  early  times  occur  while  decrease  in  wellbore  storage  is  reducing  flow  into  the  borehole  from  the  formation.  At 
later  times,  the  influence  of  wellbore  head  decline  on  flow  into  the  borehole  has  become  negligible. 

At  the  location  identified  as  EBF-2  in  Figure  5,  the  receiving  cell  of  the  observation  is  identified  as  Layer  9,  Row 
1,  and  Column  1.  The  discharging  cell  is  identified  as  Layer  10,  Row  1,  and  Column  1.  Thus,  positive  flow  is  defined 
as  upward  within  the  borehole.  A  weighting  factor  of  1.0  is  used  because  both  a  flow  meter  measurement  and  the 
internal  flow  observation  should  account  for  all  the  flow  in  the  borehole. 


Figure  6  presents  the  simulated  measurements  for  the  three  locations  identified  in  Figure  5.  The  upper  curve, 
EBF  1,  includes  all  the  flow  from  the  formation.  The  small  flows  initially  indicate  almost  the  entire  response  to 
pumping  is  a  decrease  in  the  water  level  in  the  wellbore.  As  time  progresses,  a  larger  fraction  of  the  pumping  is  made 
up  of  formation  fluid.  After  a  dimensionless  time  of  103  the  water  level  in  the  wellbore  has  stabilized  such  that  nearly 
all  of  the  pumping  rate  is  made  up  by  formation  fluid,  providing  confirmation  that  the  measurement  calculation  is 

correct. 

If  we  inspect  the  late  time  flow  rates,  we  see  that  19%  of  the  flow  into  the  wellbore  occurs  from  layer  11,  38% 
from  Layer  10,  43%  from  Layer  9.  In  the  first  test  case,  the  pumping  withdrawal  was  apportioned  between  Layers  9, 
10,  and  11  based  on  the  thickness  of  the  layers,  resulting  in  a  flow  split  of  20%,  40%,  and  40%  respectively.  Figure 
6  suggests  the  assumed  flow  split  was  in  error.  Obviously,  Figure  3  indicates  that  the  error  was  not  significant.  The 
discrepancy  highlights  a  potential  use  of  the  cylindrical  geometry  model:  determining  relative  pumping  rates  into 
layers  of  a  multilayer  well  for  use  in  larger  regional-scale  simulations. 

Column 
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Figure  5:  Detail  of  Internal  Flow  observation  locations 
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Figure  6:  Normalized  flow  at  the  internal  flow  observation  locations. 
The  locations  are  defined  in  Figure  5. 


4  Incorporation  of  Water  Table  Movement  into  a  Block- Centered  Finite- 
Difference  Model 

4.1  Introduction 

In  this  section,  we  demonstrate  that  simnlation  of  fluid  flow  near  a  water  table  is  improved  if  the  node  of  a  cell 
that  contains  the  water  table  is  located  at  the  water  table  rather  than  at  the  center  of  the  cell.  The  location  of  the 
node  should  dynamically  Mow  the  water  table.  MODFLOW-2000  uses  the  water  table  elevations  to  calculate  the 
hydraulic  conductance  but  not  to  position  the  node  location.  Since  the  water  table  elevations  are  already  available, 

a  minor  change  to  the  code  results  in  improved  performance. 

Tb  properly  develop  the  sound  rationale  for  our  modification,  we  need  to  review  some  of  the  fundamentals  of  the 
finite-difference  formulation  of  groundwater  flow  that  we  avoided  earlier  by  stressing  the  relation  of  the  cylindrical 
geometry  modification  to  the  RADMOD  code.  MODFLOW  is  based  on  a  thre^dimensional  block-centered  finite- 
difference  formulation.  Tbe  subsurface  is  divided  into  a  set  of  contiguous  cells.  The  flow  properties  of  each  cell 
are  specified  by  a  set  of  hydrologic  parameters.  The  equations  of  flow  through  the  subsurface  are  approximated 
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by  finite-difference  equations  that  use  a  single  head  value  to  describe  the  fluid  pressure  associated  with  a  cell.  The 
equations  axe  formed  using  the  assumption  that  the  head  exists  at  a  specific  location,  the  cell’s  node.  The  term 
block-centered  indicates  the  node  is  located  at  the  center  of  the  cell. 


4.2  Block- Centered  Finite-Difference  Formulation 

The  equation  for  saturated  flow  in  terms  of  hydraulic  head,  h}  is 

„  /_  „  dh(x,y,z)  _  8  _  _s dh(x, y, z)  ,  d  „  _^dh(x,y,z)  ,  d  _N dh(x,y,z )  ,  ..  /nnS 

s»(x,y,z)  ^  -  dxK(x,y,z)  ^  +  dyK(x,y,z)  ^  +  QzK(x,y,z)  &  +Q(x,y,z)  ( 20) 

where:  Sa  is  specific  storage,  K  is  hydrauhc  conductivity,  and  Q  is  a  fluid  source  or  sink.  A  finite-difference  model 
is  a  discretization  of  the  flow  equation  where  the  head  is  only  calculated  at  a  finite  number  of  locations.  Material 
properties  are  considered  constant  within  cells.  Figure  7  shows  nine  cells  of  a  block-centered  finite-difference  model. 
The  figure  represents  a  vertical  cross-section,  of  thickness  Ac  into  the  page,  with  elevations  given  by  the  z  coordinates, 
and  horizontal  location  given  by  r.  The  circles  axe  cell  nodes. 

Column 
1  2 


<5 


S  2 


°1 

*2 

°3 

*4 

#5 

*6 

AZ2 

°7 

•fc 

AZ3 

Ari 

Ar2 

Ar3 

Figure  7:  Example  block-centered  finite-difference  grid. 

Dots  and  circles  represent  node  positions. 

The  nodes  represented  as  dots  are  in  hydraulic  contact  with  the  center  node. 
The  nodes  depicted  as  open  circles  are  not. 
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To  form  the  finite-difference  equations,  the  filled  nodes  are  considered  to  be  in  hydraulic  contact  with  the  center 
cell.  The  cells  in  contact  diagonally,  open  circles,  do  not  connect  hydraulically  to  the  center  cell.  The  gradient  is 
calculated  by  treating  the  hydraulic  head  of  a  cell  to  be  located  at  the  node,  the  center  of  the  cell.  Fluid  flow  between 
two  cells  is,  in  effect,  one  dimensional.  The  one  dimensional  flow  is  converted  to  volumetric  flow  by  multiplying  by 
the  area  of  the  contact  face  between  cells.  In  Figure  7  the  contact  area  between  nodes  5  and  6  is  AcAz2.  The  flow 
into  the  face  from  node  5  is  the  same  as  the  flow  out  of  the  face  to  node  6. 


Qte  —  Q&6 

(21) 

AcA  «Kjh*-hs)  =  AcaW'"!^ 

2  2 

(22) 

CiZs  (/&56  hg)  =  CJt$  (he  —  h$e) 

(23) 

where:  Qse  represents  the  flow  across  the  contact  face  of  cells  5  and  6;  Kr  is  the  hydraulic  conductivity  in  the  r 
direction;  h5&  is  the  head  at  the  interface  of  the  cells  5  and  6;  and  CR,  using  the  MODFLOW  convention,  is  the  cell 
conductance  in  the  r  direction.  Eliminating  h66  from  (23)  we  get 

056  =  C^'^Rg  {h*  ~  ^  (24) 

=  CJ?66  (he  —  hs)  (25) 


where  CR§§  is  the  conductance  between  nodes  5  and  6  as  in  (4). 

Integrating  over  the  cell  volume,  the  left  hand  side  of  (20)  is  approximated  by 


S*AcArAz 


A  hi 
A  t 


where  S$  is  the  specific  storage  of  the  cell  material,  and  Si  is  the  capacitance  of  the  cell.  Specific  storage  represents 
the  change  in  void  space  within  the  aquifer  that  occurs  in  response  to  a  change  in  the  fluid  pressure.  As  the  head 
increases  the  load  that  must  be  supported  by  the  aquifer  material  decreases.  The  aquifer  expands  slightly  in  response 
to  decrease  in  load  creating  additional  void  space  for  fluid  storage.  The  discretized  flow  equation  for  an  internal  cell 
is 

Ah  N 

5*^T  =  E C**  fa  -  h<)  +  &  (26) 

where:  subscript  a:  is  a  meta  character  that  becomes  r,  c,  or,  v  for  flow  between  columns,  rows  or  layers  respectively; 


i 
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N  is  the  number  of  neighboring  cells,  4  for  two-dimensions  and  6  for  three  dimensions;  and  Qi  may  account  for  flow 
sources  or  sinks  to  the  node  in  addition  to  fluid  flow  through  the  aquifer. 

Having  developed  the  finite- difference  flow  equation  for  internal  nodes,  we  turn  to  a  cell  that  represents  the  top 
layer  of  an  unconfined  aquifer.  Fluid  flow  in  the  unsaturated  portion  of  the  aquifer  is  approximated  by  two  terms,  Qi 
to  represent  recharge  and  Sv  to  represent  the  effects  of  changing  water  table  elevation.  Sy  accounts  for  the  change 
in  fluid  occupied  void  space  within  the  cell  for  a  change  in  water  table  elevation.  Typically,  in  a  fully  saturated  flow 
code  such  as  MOD  FLOW,  the  dynamics  of  water  movement  within  the  capillary  fringe  are  ignored  by  treating  the 
change  in  occupied  void  space  to  be  instantaneous.  Sy  can  be  approximated  as 

Sy  =  {&  -  ©r) 

where  ©  is  porosity  and  0r  is  residual  saturation  of  the  aquifer  material.  As  the  change  in  water  table  elevation  is 
equivalent  to  the  change  in  hi ,  the  capacitance  of  a  water  table  cell  becomes 

Si  =  (S,Az  +  Sy)  AyAr 
Si  -  SyAcAr  =  SYi 

S8Az  has  been  dropped  in  (28)  because  it  is  usually  negligibly  small  compared  to  Sy. 

It  is  important  to  keep  in  mind  that  a  typical  fully  saturated  code  approximates  the  water  table  as  a  sharp 
boundary,  fully  saturated  below  the  water  table  and  completely  unsaturated  above.  Variations  to  the  hydraulic 
conductivity  above  the  water  table  based  on  proximity  to  the  water  table  (Doherty,  2001)  do  not  change  the  sharp 
boundary  approximation  for  specific  yield. 

Replacing  SSi  by  SY*  in  (26)  ignores  a  difference  between  specific  storage  and  specific  yield.  I  find  a  conceptual 
sleight  of  hand  helps  illustrate  the  difference.  The  concept  of  a  physical  change  in  the  void  space  is  replaced  by  the 
concept  of  fluid  creation  or  annihilation.  The  concept  is  consistent  with  the  phrase  ‘Svater  coming  out  of  storage” 
that  is  often  used  to  describe  the  effects  of  specific  storage.  The  equations  of  flow  would  be  identical  if  instead  of  a 
decrease  in  void  space  an  equivalent  increase  in  the  volume  of  fluid  occurred.  With  respect  to  specific  storage,  this 
source  of  creation  and  removal  of  water  would  be  distributed  throughout  the  cell  volume.  To  approximate  the  effect 
with  a  node  at  the  center  of  the  cell  is  appropriate. 

In  the  case  of  specific  yield,  the  change  in  available  void  space  to  the  fluid  all  occurs  above  the  water  table, 
approximated  in  a  fully  saturated  flow  code  as  a  distinct  boundary.  For  a  change  in  the  water  table  position,  the 


(27) 

(28) 
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fluid  source  should  be  positioned  at  the  water  table  for  it  is  there  that  the  change  in  void  space  occurs.  Approximating 
the  source  at  the  center  of  the  cell  introduces  an  unnecessary  error. 

We  now  look  at  the  impact  of  positioning  a  node  at  the  water  table  on  the  equations  of  flow  for  cells  that  include 
the  water  table.  Figure  8  depicts  three  cells.  The  water  table  is  represented  by  the  dashed  line  in  the  upper  cell. 
Large  dots  are  placed  at  the  center  of  each  cell.  A  circle  is  located  at  the  center  of  the  water  table  within  the  upper 
cell.  For  a  change  in  the  head  of  the  middle  cell,  the  head  dynamics  of  the  top  cell  are  dominated  by  the  flow  across 
the  boundary  between  the  two  cells  balanced  by  the  creation  of  water  at  the  water  table.  Location  of  the  cell  node 
is  incorporated  in  the  vertical  conductance  of  the  upper  cell.  If  the  node  is  located  at  the  center  of  the  cell,  then 


ck 


AcArKvi 


Azi 

2 


(29) 


The  superscript  c  refers  to  the  node  at  the  center  of  the  cell.  Kv\  is  the  vertical  hydraulic  conductivity  of  cell  1.  A 
superscript  wt  is  used  to  indicate  the  node  at  the  water  table.  If  the  node  is  located  at  the  water  table  then 


AcArKvi 
zwt  —  (zz  +  Af*) 


(30) 


or 


AcArKvi 
hi-(z»  +  *p) 


T 

AZ. 

+ 

AZ2 

+ 

AZ, 

JL 


(31) 


Figure  8:  Cells  with  water  table. 
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If  the  water  table  exists  at  the  top  of  the  cell  then  (31)  becomes 


Cwt 

vl 


AcArKvi 

Azi 

Cl , 

2 


(32) 

(33) 


Likewise,  the  conductance  between  the  upper  and  middle  cells  can  be  related  as 


Cwt 
v!2 


1  Cgi  +  Cv  2 
+  Cv  2 


2Si 

2 


^12 


(34) 


For  the  case  of  Azi  —  Az2  and  Kv i  —  KV2> 


Cwt  ^  /tc 

vl2  ““  g°t>12 


(35) 


For  Figure  8,  if  the  water  table  is  near  the  bottom  of  the  cell  rather  than  near  the  top,  the  conductance  between 
nodes  1  and  2  becomes 


r*wt  _ one 

°vX2  — 


v!2 


(36) 


4.3  Implementation 

Only  a  small  change  in  coding  is  needed  to  implement  a  water-table-following  node.  In  the  subroutine  SGF1LPF1VCOND 
(gwfllpfl.f),  the  variable  BOVK1  represents  £  for  the  upper  cell  which  we  can  relate  to  (29)  by  C&  =  .  The 

calculation  of  BOVK1  changes  from 

B0VK1  =  0.6*  (B0TM(  (  J ,  I  ,LB0TM(K)  - 1)  -B0TH(  ( J , I ,LB0TM(K)  )  )  /HYC1 

for  a  normal  cell,  to 


B0VK1  =  (HNEW(I , J ,K) -B0TM( ( J , I , LBOTM(K) ) ) /HYC1 

when  the  head  (HNEW)  in  the  cell  is  less  than  the  top  of  the  cell  (BOTM((J,I,LBOTM(K)-l)).  BOTM((J,I,LBOTM(K))) 
is  the  bottom  of  the  cell  and  HYC1  is  the  vertical  hydraulic  conductivity,  Kv,  of  the  cell. 
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4.4 


Test  Cases 

4.4.1  Three-Dimensional  Comparison  to  WTAQ 

The  influence  of  the  position  of  the  water  table  node  was  compared  to  WTAQ  calculations  for  a  hypothetical 
homogeneous  aquifer.  The  geometry  and  properties  of  the  aquifer  and  the  positions  of  the  monitoring  wells  are 
similar  to  the  Boise  Hydrogeophysical  Research  Site,  but  the  similarity  is  not  significant  to  this  study.  The  aquifer 
properties  are  listed  in  Table  1. 


Parameter 

Value 

Horizontal  K 

l.xl0*4  m/s 

Vertical  K 

l.xl0'4  m/s 

Specific  storage 

l.xlO*4  per  m 

Specific  yield 

0.38 

Aquifer  thickness 

18  m 

Table  1:  Aquifer  characteristics  for  the  three-dimensional  water  table  modification  test. 

The  pumping  well  has  a  radius  of  0.0508  m  and  is  screened  over  an  interval  of  8  m  to  12  m  below  the  water  table. 
Two  MODFLOW  models  using  Cartesian  geometry  are  tested.  Both  have  the  same  43  rows  and  45  columns.  The 
horizontal  grid  spacing  was  constructed  such  that  nodes  occurred  dose  to  the  following  distances  from  the  pumping 
well:  3.52  m  ,  7.08  m,  and  22.1  m.  One  model  has  21  layers  the  other  26.  Each  layer  in  the  21  layer  model  is  1  m 
thick  except  for  0.5  m  layers  above  and  below  each  end  of  the  pumping  interval,  and  at  the  top  and  bottom  of  the 
model.  The  26  layer  model  has  a  fine  grid  near  the  water  table.  In  the  26  layer  model,  the  uppermost  layer  is  0.1  m 

thick  with  eight  layers  in  the  upper  3  meters. 

Figure  9  presents  calculations  of  drawdown  at  distances  from  the  pumping  well  of  3.52  m  (well  B2),  7.08  m  (well 
C3),  and  22.1  m  (well  X4)  at  depths  of  1  m,  6  m,  and  10  m  below  the  water  table.  Four  sets  of  data  are  presented 
in  each  panel.  All  of  the  MODFLOW  simulations  closely  match  the  WTAQ  calculations.  The  largest  deviation 
between  the  calculations  occurs  at  nodes  closest  to  the  water  table  and  at  times  when  the  drawdown  is  first  impacted 
by  the  water  table  decline,  causing  a  brief  plateau  in  the  drawdown.  After  the  plateau  region  the  impact  of  water 
table  decline  is  fast  with  respect  to  the  time  scale.  Figure  10  shows  the  region  of  greatest  disagreement,  the  plateau 
region  of  the  location  nearest  to  the  water  table  and  closest  to  the  well.  In  this  figure,  the  differences  between  model 
calculations  are  clearly  presented.  The  dot-dash  line  of  the  26  layer  model  is  a  definite  improvement  over  the  dashed 
line  of  the  21  layer  model.  The  dotted  line  of  the  modified  MODFLOW  code  fails  to  overlay  the  WTAQ  calculations 

only  at  very  early  times. 
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Figure  9:  Comparison  of  drawdown  at  various  locations. 

WTAQ  calculation  (solid),  unmodified  MODFLOW  with  21  layers  (dashed),  unmodified  MODFLOW  with  26  layers 
(dot-dash),  modified  MODFLOW  with  21  layers(dotted).  The  dotted  line  overlays  the  solid  line  at  these  scales. 
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For  the  test  case,  the  water  table  modification  slows  the  execution  speed  of  the  model  by  7%.  The  increase  in 
layers  to  26  compared  to  21  layers  requires  31%  more  time  to  reach  a  solution. 


Location  B2-1m 


Figure  10:  Comparison  of  drawdown  at  a  depth  of  1  m  below  the  water  table  3.52  m  away  from  the  well. 

WTAQ  (solid),  unmodified  with  21  layers  (dashed),  unmodified  with  26  layers  (dot-dash),  and  modified  with  21 

layers  (dotted). 


4.4.2  Two-Dimensional  Cylindrical  Geometry  Comparison  to  WTAQ 

We  use  the  test  case  of  Section  2.3  for  the  cylindrical  geometry  modification.  Figure  11  presents  the  results  of  the 
cylindrical  geometry  test  case  with  the  modified  water  table  simulation.  WTAQ  calculations  are  presented  as  a  solid 
line.  The  cylindrical  model  calculations,  using  the  center  of  the  uppermost  layer  as  the  node  position,  Eire  shown 
with  a  dashed  line.  The  dotted  line  presents  the  calculations  with  the  uppermost  node  following  the  water  table. 
The  three  node  positions  are  located  by  layer  and  radial  column  number.  For  example,  Node  3  2  is  positioned  in 
Layer  3  Column  2  which  translates  to  20  ft  below  the  water  table  and  0.52  ft  away  from  the  well. 
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Figure  11:  Comparison  of  WTAQ  calculations  (solid),  cylindrical  model  without  water  table  modification  (dashed), 
and  modified  cylindrical  model  for  the  RADMOD  test  case. 

5  Limitations  of  code  modifications 

This  section  details  the  known  limitations,  and  obvious  potential  limitations,  in  the  current  status  of  our  modification 
to  MODFLOW-2000.  The  limitations  are  with  respect  to  compatibility  with  all  other  packages  of  MODFLOW-2000. 
The  compatibility  to  most  other  packages  have  yet  to  be  investigated.  Incompatibilities  that  have  been  identified 
can  be  rectified  if  need  be.  The  need  for  a  FORTRAN  90  or  95  compiler  could  be  removed,  but  I  prefer  the  present 
organization  of  the  code. 

The  layer  property  flow  package  has  been  modified  to  accommodate  cylindrical  geometry,  but  the  block-centered 
flow  package  has  not  been  modified.  The  implication  of  the  lack  of  block-centered  flow  package  compatibility  is  that 
the  LPF  input  file  must  be  used  rather  that  the  BCF6  input  file  to  describe  the  material  properties  in  the  model. 
The  horizontal  flow  barrier  package  has  been  modified  to  accommodate  cylindrical  geometry,  but  the  river,  recharge, 
drain,  evapotranspiration,  ADV,  and  general-head  boundary  packages  have  not  been  investigated  with  respect  to 
compatibility.  It  is  anticipated  that  none  of  these  packages  will  be  difficult  to  update  for  the  cylindrical  geometry. 
It  may  be  that  most  require  no  changes. 

We  have  not  investigated  the  impact  of  the  water  table  modification  on  the  sensitivity  process.  There  clearly  is 
an  impact  that  will  most  likely  need  to  be  addressed  in  a  manner  similar  to  the  influence  of  water  table  variation  on 
horizontal  conductance.  We  intend  to  address  this  limitation  soon. 
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6  Summary 

This  report  documents  three  separate  enhancements  to  the  MODFLOW-2000  code:  cylindrical  geometry,  interna) 
flow  measurements,  and  improved  water  table  simulation.  The  cylindrical  geometry  enhancement  allows  the  use  of 
the  observation,  sensitivity  and  parameter  estimation  capabilities  of  MODFLOW-2000  to  be  applied  to  situations 
where  radial  symmetry  can  be  assumed.  Potential  uses  of  cylindrical  geometry  are  the  simulation  of  aquifer  dynamics, 
estimation  of  the  influences  of  disturbed  zones  around  a  wellbore,  determination  of  flow  splits  into  layers  that  can 
be  applied  to  larger  scale  simulations,  and  our  intended  use  -  the  interpretation  of  borehole  flow  meter  data. 

The  addition  of  internal  flow  observations  to  the  observation,  sensitivity  and  parameter  estimation  processes  is  a 
relatively  minor  enhancement.  In  addition  to  simulation  of  flows  for  parameter  estimation,  the  simulated  observations 
could  useful  as  an  alternative  to  obtaining  flow  rates  from  the  cell-by-cell  flow  file.  The  enhancement  is  needed  for 
our  intended  investigation  of  borehole  flowmeter  interpretation. 

Our  change  in  water  table  simulation  will  improve  the  performance  of  any  MODFLOW  simulation  when  water 
table  dynamics  are  important.  This  improvement  has  been  explained  in  terms  of  the  process  of  water  table  rise  and 
drainage.  In  the  limited  test  cases  we  have  run,  the  modifications  lead  to  instability  in  the  model  equations  that 
resulted  in  less  than  a  10%  increase  in  computation  time  and  in  noticeable  improvement  in  accuracy. 
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A  Appendix  A 

Input  Instructions 

A.l  Input  instructions  for  inner  boundary  of  radial  geometry 

The  file  contains  a  single  number,  the  distance  of  the  inner  boundary  of  the  radial  cross-section  from  the  axis  of 
symmetry.  The  file  name  is  specified  in  the  name  file  using  file  type  designation  CGEO. 

0.  [text] 

Item  0  is  optional  and  can  include  as  many  lines  as  needed. 

Each  line  must  contain  the  #  character  in  first  column. 

1.  SRI  (free  format) 


Explanation  of  Variables 

text  is  a  character  string  of  up  to  79  characters  that  starts  in  column  2. 

SRI  is  the  distance  from  the  axis  of  symmetry  in  the  cylindrical  geometry  to  the  inner  boundary  of  Column 

1.  The  column  widths  are  defined  by  DELR  of  the  DIS  file.  The  columns  progress  outward  from  SRI. 

A. 2  Input  instructions  for  internal  flow  observations 

Input  instructions  for  the  internal  flow  package  are  read  from  a  file  specified  with  FLOB  as  the  file  type  in  the  name 
file. 

0.  [#text] 

Item  0  is  optional  and  can  include  as  many  lines  as  needed.  Each  line  contains  the  #  character  in  first 
column. 

1.  NQFL  NQCFL  NQTFL  (free  format) 

2.  TOMULTFL  EVFFL  IOWTQFL  (free  format) 

Read  items  NQTFL, 4,  and  5  for  each  NQFL  groups  of  cells  for  which  internal  flow  observations  are  to  be 
specified. 

3.  NQOBFL  NQCLFL  (free  format) 

Read  item  4  for  each  NQOBFL  observation  times  for  this  group  of  cells.  STATISTIC  and  STAT-FLAG 
are  ignored  if  IOWTQFL  is  greater  than  zero. 
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4.  OBSNAM  IREFSP  TOFFSET  HOBS  STATSTIC  STAT-FLAG  PLOT-SYMBOL  (free  format) 

Read  item  5  for  each  of  |NQCLFL|cells  in  this  group. 

5.  TO _ LAYER  TO_ROW  TO_COL  FR_LAYER  FR_ROW  FR_COL  FACTOR  (free  format) 

Read  items  6  and  7  if  IOWTQFL  is  greater  than  zero. 

6.  FMTIN  IPRN  (free  format) 

7.  WTQ(1,1),  WTQ(1,2),  WTQ(1,NQTFL),  ....  WTQ(1,NQTFL)  (format:  FMTIN) 

WTQ(2,1),  WTQ(2,2),  WTQ(2, NQTFL), WTQ(2,NQTFL) 


WTQ(NQTFL,1),  WTQ(NQTFL,2),  WTQ(NQTFL,3),  WTQ(NQTFL, NQTFL) 


Explanation  of  Variables 

text  is  a  character  string  of  up  to  79  characters  that  starts  in  column  2. 

NQFL  is  the  number  of  cell  groups  for  which  internal  flow  observations  are  listed.  A  group  in  a  collection  of  cell 

pairs  that  define  an  internal  surface  that  represents  the  location  of  a  single  flow  measurement. 

NQCFL  is  greater  than  or  equal  to  the  total  number  of  cells  in  all  cell  groups.  NQCFL  must  be  greater  than  or 
equal  to  the  sum  of  all  |NQCLFL| 

NQTFL  is  the  total  number  of  internal  flow  observations  for  all  cell  groups.  NQTFL  must  equal  the  sum  of  all 
NQOBFL,  which  are  specified  in  repetitions  of  item  3. 

TOMULTFL  is  the  time  offset  multiplier  for  the  observations.  The  product  of  TOMULTFL  and  TOFFSET  must 
produce  a  time  value  in  units  consistent  with  other  model  input.  TOMULTFL  can  be  dimensionless  or 
can  be  used  to  convert  the  TOFFSET  units. 

EVFFL  is  the  error  variance  multiplier  for  internal  flow  observations.  It  is  used  to  calculate  the  weights  as 
described  in  the  STATSTIC  explanation. 

IOWTQFL  indicates  whether  the  variance-covariance  matrix  is  to^be  read  for  the  internal  flow  observations  as  item 
7.  If  IOWTQFL  is  greater  than  zero,  then  the  variance-covariance  matrix  is  read  and  used  to  calculate 
the  weights.  If  IOWTQFL  is  zero,  the  matrix  is  not  read  and  STATISTIC  is  used. 
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NQOBFL  is  the  number  of  times  at  which  flows  are  observed  for  the  group  of  cells. 


NQCLFL  is  used  for  two  purposes.  The  absolute  value  of  NQCLFL  is  the  number  of  cell  pairs  in  the  group.  If 
NQCLFL  is  negative  then  FACTOR  is  forced  to  1.0  for  all  cells  in  the  group,  otherwise  FACTOR  is  read 
for  each  cell  pair  in  item  5. 

OBSNAM  is  a  1  to  12  character  string  used  to  identify  the  observation.  OBSNAM  may  not  contain  internal  white- 
space  characters. 

IREFSP  is  the  reference  stress  period  to  which  the  observation  time  is  referenced.  The  reference  time  is  the 
beginning  of  the  stress  period. 

TOFFSET  is  the  time  from  the  beginning  of  the  stress  period  IREFSP  to  the  time  of  the  observation.  The  multiplier 
TOMULTFL  allows  TOFFSET  to  be  in  units  other  than  the  units  of  the  simulation. 


HOBS  is  the  observed  volumetric  flow  rate  across  the  boundary  defined  by  the  group  of  to  and  from  cell  pairs. 

STATISTIC  is  the  value  from  which  the  weight  for  the  observation  is  calculated.  STAT-FLAG  controls  the  calculation 
as  described  below.  STATIC  is  ignored  if  the  variance-covariance  matrix  is  used  (IOWTQFL>0). 


STAT-FLAG  controls  how  STATICSTIC  is  used.  If  IOWTQFL>0  then  STAT-FLAG  is  ignored 


STAT-FLAG=0  =>  STATISTIC  is  a  scaled  standard  variance 


(  L*\‘ 
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weight  =  (STATl  STIC.  £v  FPL) 

STAT-FLAG=1  =>  STATISTIC  is  a  scaled  standard  deviation 
weight  =  (sTATlSri&.BVPPrj 

STAT-FLAG=2  =>  STATISTIC  is  a  scaled  coefficient  of  variation  [] 
weight  =  (statistic*hobs)2*evffl 

PLOT-SYMBOL  is  an  integer  that  will  be  written  to  the  output  files  intended  for  graphical  analysis  to  allow  control 
of  the  symbols  used  when  plotting  data. 


TO-LAYER  is  the  layer  of  the  recieving  cell.  Positive  flow  is  toward  the  recieving  cell  and  away  from  the  discharging 
cell.  TO-LAYER  and  FR-LAYER  define  the  layer  boundary  across  which  the  simulated  flow  measurement 
is  calculated. 


TO-ROW  is  the  row  of  the  recieving  cell.  TO-ROW  and  FR-ROW  define  the  row  boundary  used  to  calculate  the 


simulated  flow  measurement. 


TO-COL  is  the  column  of  the  recieving  cell.  TO-COL  and  FR-COL  define  the  column  boundary  used  to  calculate 
the  simulated  flow  measurement. 


FR-LAYER  is  the  layer  of  the  discharging  cell. 

FR-ROW  is  the  row  of  the  discharging  cell. 

FR-COL  is  the  column  of  the  discharging  cell. 

FACTOR  is  the  fraction  of  the  calculated  flow  added  to  the  flow  measurement. 

FMTIN  is  the  format  used  to  space  the  variance-covariance  matrix  data  listed  in  item  7.  The  format  is  a  FOR¬ 
TRAN  specification  for  REAL  numbers  and  must  be  enclosed  in  parentheses. 


IPRN  identifies  the  format  used  in  printing  the  variance-covariance  matrix.  If  IPRN  is  less  than  1  the  matrix 
is  not  printed. 


Output  requires  more  than  80  columns 

IPRN 

FORMAT 

1 

10G12.3 

2 

10G12.4 

3 

9G12.5 

4 

8G13.6 

5 

8G13.7 

Output  fits  in  80  columns 

IPRN 

FORMAT 

6 

5G12.3 

7 

5G12.4 

8 

5G12.5 

9 

4G13.6 

10 

4G13.7 

Table  2:  IPRN  format  specifications 


WTQ  Is  an  NQTFL  by  NQTFL  array  containing  the  variance-covariance  array  for  the  internal  flow 
measurements  in  units  of  [?j%l .  All  elements  of  the  matrix  must  be  entered. 
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